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I. INTRODUCTION 

Biological molecules self-assemble into membranes, protein assemblies, viruses and cells.- 
Material design inspired by nature is a promising route to create materials with novel or 
enhanced properties by spontaneous self-assembly.-'- In the laboratory, colloidal particles 
can be synthesized with a variety of shapes and directional interactions.- These patchy 
particles could potentially be used to mimic the self-assembly observed at smaller length 
scales, and to rationally design assemblies from their basic building blocks.- 

Studies of self-assembly range from those considering only repulsive interactions which 
dehne the shape of the particle,->^ to those considering spherical particles with directional 
attractions.—'— In colloidal systems, both the shape and the directional interactions are 
intimately coupled when depletant is added to the solution.— This depletant interaction 
drives the assembly of lock-and-key colloids.— For colloidal clusters synthesized with smooth 
and rough beads, the smooth beads attract more strongly to one another than to rough 
beads, due to more excluded volume overlap at contact.— The focus of this paper is on the 
self-assembly of trimers, consisting of a central attractive bead and two repulsive beads. 

In an experimental and computational study, it was observed that dimers with one at¬ 
tractive bead and one repulsive bead self-assembled into spherical micelles.— In addition, 
it was demonstrated that trimers with one attractive bead and two repulsive beads could 
be synthesized. Recently, an experimental and computational study of trimers with one 
attractive bead and two repulsive beads reported that only elongated clusters were formed, 
in contrast to the spherical micelles formed by dimers.— In a different computational study, 
a flexible 3-mer chain with two attractive beads and one repulsive bead at the end of the 
chain was found to self-assemble.— Tetramers with two attractive beads and two repulsive 
beads were found to self-assemble into a variety of structures and used to study protein 
aggregation.— 

In this computational study, the phase behavior of a family of trimer models with one 
attractive central bead and two repulsive beads is investigated for a range of different trimer 
shapes. Advanced simulation methods were used to obtain the fluid phase behavior based 
upon thermodynamic and structural dehnitions, rather than more phenomenological ap¬ 
proaches. In particular, Wang-Landau Transition Matrix Monte Carlo (WL-TMMC) simu¬ 
lations were preformed in the grand canonical ensemble, utilizing on the order of hundreds 
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of billions of trials per simulation. The trimers form spherical micelle-like clusters, elongated 
clusters and densely packed cylinders. We show that there is a transition from self-assembly 
to bulk fluid phase separation as bond length is reduced, and hnd the in-between bond 
length with both elongated, self-assembled structures and fluid phase separation. We also 
discuss how the phase behavior of the family of trimer models may be understood in terms 
of the interaction between particles. 

This paper is organized as follows. In Sec. [TTl we describe the family of trimer models 
studied in this work. We then discuss the computational methods and the thermodynamic 
and structural transition dehnitions in Sec. m Results are discussed in Sec. HVl Finally, 
we conclude and discuss future work in Sec. El 


II. MODELS 


In this paper, we studied the fluid phase behavior of a family of trimer models. The trimer 
consisted of one central, attractive bead and two repulsive beads, as shown in Figure [Hand 
Table [B Specihcally, we studied how fluid phase behavior was affected by simple geometric 
parameters, L and 0, where L is the rigid bond length between an attractive and repulsive 
bead, and 0 is the rigid bond angle with the vertex on the attractive bead. The beads 
interact via a shifted-force Lennard-Jones (LJ) potential. 


UlJ(r) = ULj{r) - ULj{rc) - (r - rj 


dU, 


LJ 


dr 


ULj{r) = 46 


a 


12 


r=r^ 

CT\ 6 


( 1 ) 

( 2 ) 


where Tc is the potential cut-off, U{r > Vc) = 0. For interactions between attractive, blue 
beads, rdcr = 3. All other pair-wise interactions are purely repulsive, rdo = 2^/®, also 
known as the Weeks-Chandler-Andersen potential.— Each bead has equal a, e and mass. 

While the model described above may seem simplistic, it is intended to capture basic 
geometric features that should be relevant to a broad range of systems. Indeed, this trimer 
model exhibited rich phase behavior with respect to self-assembly and fluid phase separation 
(see Figures [2] and |3]). The aim of this study is to rationalize how the phase behavior and self- 
assembly changes with particle shape, using a general model that may be applied to many 
different types of systems and is computationally tractable. In this study, the Lennard-Jones 
potential was chosen for simplicity. 
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FIG. 1. The family of trimer models investigated in this work illustrated using VMD.— Blue 
beads are attracted to other blue beads, while all other pair interactions (red-red and blue-red) 
are purely repulsive. The trimers are listed in order of increasing attractive region with respect to 
repulsive region, and the same order as Table HI 


TABLE I. Trimer model parameters, L and 0, and computed values for the excluded volume (see 
Appendix [B]), critical temperature, and Boyle temperature (see Appendix IG]). 


Lja 

0 

CO 

ksTc/e 

fesTeoyie/ e 

1 

7r/2 

9.83 

n/a 

0.365(5) 

1 

7r/3 

9.31 

n/a 

0.435(5) 

1 

7r/4 

8.88 

n/a 

0.485(5) 

0.75 

7r/3 

8.02 

n/a 

0.505(5) 

0.4 

7r/3 

6.19 

0.3117(1) 

0.815(5) 

0.25 

7r/3 

5.41 

0.4989(1) 

1.17(1) 

0 

7r/3 

4.19 

0.8798(7) 

2.00(2) 


III. METHODS 

Flat-histogram sampling methods were used to investigate the fluid phase behavior of the 
family of trimer models. Specihcally, Wang-Landau Transition-Matrix Monte Carlo (WL- 
TMMC) simulations^^— in the grand canonical ensemble were performed, as described 
below in Section IIII A1 This powerful simulation method computes the free energy, potential 
energy and pressure as a function of density at constant temperature, as well as provides 
detailed structural information, in a single simulation. The advantage of the grand-canonical 
ensemble over the canonical ensemble is that smaller system volumes can be used to capture 
physically relevant density fluctuations. In the canonical ensemble, where the total num¬ 
ber of particles is hxed, the use of small system volumes amounts to the imposition of a 
constraint.— For self-assembling systems, this arbitrary constraint may not agree with the 
preferred free monomer densities and sizes of self-assembled structures in the thermodynamic 
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FIG. 2. Illustration of selected structures. Unless otherwise specified, L = a, Q = tt/S, ksT/e = 
0.2, V = 729(T^ and the blue boxes represent periodic boundaries, (a) micelle, A = 13 (b) large 
micelle, 0 = vr/d, A = 20 (c) small micelle, 0 = tt/2, A = 8 (d) elongated cluster, ksT/e = 0.125, 
A = 59 (e) elongated cluster, L = 0.4cr, ksT/e = 0.15, A = 112 (f) elongated cluster in liquid, 
L = OAa, ksT/e = 0.15, pVex = 2.6 (g) cylinder, pVex = 2.2 (h) same as (g) with a top-down 
projection 


limit. The effect of this constraint in the canonical ensemble diminishes with system size, 
and thus appropriate canonical-ensemble simulations of self-assembly require signihcantly 
larger systems that are computationally more expensive. In addition, to improve sampling 
at low temperature, WL-TMMC was combined with parallel-tempering. A single isotherm 
simulation was typically composed of hundreds of billions of trials. Simulation details are 
provided in Section IIII A1 and the methods used to determine phase coexistence and locate 
structural transitions are described in Section IIII Bl 
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A. Grand Canonical Wang-Landau Transition-Matrix Monte Carlo 


The Grand canonical WL-TMMC method was used to obtain the macrostate probability 
distribution, n(A;/i, 1/,T), which is the probability to observe the number of trimers, N, 
for a given chemical potential, /i, volume, V, and temperature, T. See Appendix A of 


Ref. 


25l for implementation details of WL-TMMC used here. The Wang-Landau update 


factor was initially set to unity, and was multiplied by 0.5 whenever the flatness criteria 
of 80 % was met. After the update factor was smaller than 10“®, the collection matrix 
was updated. After the update factor was smaller than 5 x 10“®, transition-matrix Monte 
Carlo was performed with an update to the biasing function every 10® trials. To parallelize 
the single isotherm simulations, each isotherm was divided into 12 overlapping windows to 
span the entire density range, that is the range of N values from 0 to Nmax- Because lower 
density simulations are faster, the parallelization was load balanced by decreasing window 
size with increasing density by a power scaling with exponent of 1.5. Neighboring windows 
overlapped by up to three trimers. The free energy of the entire density range was recovered 
by setting the free energy of neighboring windows equal at the middle overlapping trimer 
number, and discarding the largest and smallest number of trimers (when neighbor present) 
in each window. A window was converged if it swept at least 10 times, although most 
windows swept 100-1000 times while waiting for the high density window to converge. A 
sweep was dehned as satisfying the condition that each macrostate had been visited from 
a different macrostate at least 100 times. After a simulation swept more than one time, 
canonical ensemble averages, as described below, were accumulated for quantities such as 
the potential energy and the squared potential energy upon every successful trial attempt. 
Ideally, when running a WL-TMMC simulation, a value of the chemical potential is chosen 
such that the difference between neighboring macrostates in the macrostate distribution is 
minimized. Fortunately, the exact choice of the /i in the WL-TMMC simulations is relatively 
unimportant, because the initial Wang-Landau part of the simulation efficiently hnds the 
order of magnitude of the macrostate distribution, and then the macrostate distribution 
may be histogram reweighted to different values of /r. 


The following Monte Carlo trials were employed, as listed in Table UTl Rigid trimer transla¬ 
tions or rotations about the center of mass were attempted with equal probability. Random 
insertions or deletions of trimers were also attempted, subject to Metropolis acceptance 
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criteria.— Collective trial moves were also implemented to facilitate convergence with self- 
assembled structures. Smart Monte Carlo was used to bias the movement of trimers in 
the direction of their center-of-mass forces.— A second collective move type entailed rigid 
translation or rotation of each cluster of trimers. Clusters were dehned as all trimers having 
an attractive bead within a cut-off distance, 4 (t/ 3, from at least one other attractive bead 
in the cluster, obtained via recursive flood-hll algorithm. To obey detailed balance, cluster 
moves which resulted in a trimer joining a different cluster were rejected. Statistics on the 
clusters were accumulated every attempted cluster move, after the simulation swept more 
than one time. For each Monte Carlo trial that involved movement of trimers, the parame¬ 
ter associated with the maximum possible translation or rotation was optimized, via a 5 % 
change every 10® trials, to yield approximately 25 % acceptance of the trial move. Another 
Monte Carlo trial involved conhguration swaps between neighboring density windows to fa¬ 
cilitate convergence. Conhgurational swap moves between adjacent density windows were 
used to ensure self-assembled structures were sampled in multiple windows. These conhgu¬ 
rational swap moves helped to improve the parallelization efficiency, and were performed at 
hxed N, V, T, and /i.— Even with the assortment of trial moves described above, structural 
transitions between different self-assembled motifs were difficult to sample at low T. To 
circumvent this difficulty, parallel tempering was implemented to swap conhgurations be¬ 
tween neighboring temperatures, at hxed number of trimers, from a series of closely spaced 
isotherms. The second type of conhgurational swap move improved sampling of structural 
transitions that occurred as temperature is decreased, and was performed at hxed N, V with 
varying U, T, fi. When the conhguration swap trial is attempted, there is a 50 % chance to 
store the current conhguration, and a 50 % chance to swap the current conhguration with a 
stored conhguration on an overlapping processor (if exists), subject to Metropolis acceptance 
criteria.— 

Grand canonical ensemble averages of an observable. A, denoted as {A)^vt, are obtained 
as a continuous function of {N)^vt- Calculation of {A)^vt is based on the canonical av¬ 
erage of property A, denoted as (A) nvt, which can be calculated during the course of the 
simulation. 


(A(7V)) 


- N) 


NVT 


where rij is the number of trimers in trial i, 6 is the delta function, and N^riai are the number 


( 3 ) 
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TABLE II. Monte Carlo trials and weights for the probability of selection. 


trial 

weight 

single-trimer translation or rotation 

I 

single-trimer insertion or deletion 

1/4 

smart Monte Carlo2I 

l/WN^ax 

cluster translation or rotation 

\/hNjYiaX 

parallel configuration swap 

5 X 10"® 


of sampled states in the simulation. It follows that the grand-canonical average is 

^max 

{A{{N)f,vT))^,VT = {A{n)) NVTA{n; fi). (4) 

The quantity (AI)^v't can be obtained directly from the macrostate distribution as 

Nmax 

{N)f,vT = nn(n;/i). (5) 

n=0 

The average properties at other state conditions, namely different values of /i, can be ob¬ 
tained via histogram reweighting the macrostate distribution. 

B. Determining Phase Coexistence and Structural Transitions 

In this section, we discuss the methods used to determine fluid phase behavior. The 
two distinct types of behavior observed in this work are macroscopic phase separation and 
self-assembly (e.g. micellization), which also includes transitions between different struc¬ 
tures. Note that structural transitions that take place on a microscopic length scale, such as 
micellization, are not true thermodynamic phase transitions.— Phase coexistence conditions 
between two macroscopic phases were obtained by histogram reweighting the macrostate dis¬ 
tribution to a value of fi such that the probabilities of observing each phase are equal. Critical 
points were obtained by htting saturation densities to the law of rectilinear diameters.— 
For the remainder of this section, we discuss the methods used to determine the struc¬ 
tural transitions involving the spherical micellar fluid, which requires locating the low and 
high density and low and high temperature boundaries for the spherical micellar fluid. The 
low density boundary is dehned primarily by the critical micelle concentration, which is 
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the concentration above which a free trimer fluid becomes a micellar fluid. Similarly, the 
high temperature boundary is taken to be the maximum temperature at which micelles are 
stable. At low temperature, the micellar fluid transforms from roughly spherical structures 
to elongated ones. The temperature at which this occurs is taken to be the low-temperature 
boundary of the spherical micellar fluid. The high density transition of micelles into more 
solid-like structures is also approximately obtained, although proper sampling at high densi¬ 
ties is beyond the scope of this work. Examples of these transitions are provided in Appendix 

E 

The critical micelle concentration (CMC), dehned as the lowest concentration at which 
micelles can form, was obtained by both thermodynamic and structural dehnitions. The 
structural method directly measures the CMC as the concentration of free trimers and 
premicellar aggregates as a function of density.— This direct measurement of the CMC is 
possible because the concentration of free trimers and premicellar aggregates, pfreei remains 
approximately constant as the fluid density, {N)/V, increases at hxed temperature after 
micelle formation.^® Premicellar aggregates are dehned as clusters with a number of trimers 
less than or equal to the hrst minimum in the histogram of aggregate size (typically 4-5 
trimers), and clusters are dehned in Section UlI A1 {pfree)iiVT is relatively constant over a 
huid density range. In practice, the density range over which pfree is constant was dehned 
as the range where Pfree was within some tolerance of the hrst local maximum. In this work, 
we used a 75 % tolerance. The structurally based critical micelle concentration was taken as 
the average Pfree in fhis huid density range, and the high density boundary of the micellar 
huid was taken as the maximum density in this huid density range. The thermodynamic 
method to obtain the CMC uses the density at which the equation of state deviates from 
ideal trimer huid behavior.— The deviation appears as a second linear regime, due to the 
formation of micelles, and the density where the deviation occurs is dehned by the point 
of intersection of hts to the linear regimes. The equation of state is obtained from the 
macrostate distribution, by reweighting it to various p values, and computing the pressure 
as a function of {N)hvt/V by comparing the probability to observe zero trimers to the ideal 
trimer huid state.— In order to precisely obtain the equation of state in the density range 
of interest with WL-TMMC simulations in the grand canonical ensemble, both V and N^ax 
must be tuned to sample both the ideal and micellar huids. Depending on the temperature, 
V/a^ ranged from 9® to 64® while N^ax was 50 to 150. These low density simulations for the 
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thermodynamic dehnition of the CMC were separate from the higher density simulations 
which were used to obtain the CMC by the structural dehnition. 

The critical micelle temperature (CMT) was taken to be the highest temperature at which 
micelles could exist. This temperature is not a true critical point, and was simply named by 
analogy to the critical micelle concentration.— As noted in previous work, dehning the CMT 
is somewhat arbitrary.— Although one may dehne the CMT with structural information, it is 
difficult to distinguish between self-assembled micelles and supercritical clusters, similar to 
those formed in typical homogeneous huids. In previous work, a thermodynamic signature 
of micellization was a system-size dependent density of a second peak in the macrostate 
probability distribution of the number of particles.— Physically, the second peak corresponds 
to the formation of a micelle, which happens at the same number of trimers, regardless if 
the system size is slightly increased. This thermodynamic signature of micellization was 
used to dehne the CMT in this work, as demonstrated in Appendix |Al The error bars for 
the CMT simply depended on the spacing between simulated isotherms. Specihcally, the 
CMT was obtained by identifying two isotherms. The hrst is the highest temperature in 
which the huid contained micelles, and the second is the next highest temperature in which 
the huid did not contain micelles. This ehectively brackets the CMT. Therefore, the CMT 
must be in between these two isotherms. The reported value of the CMT was the average 
of these two isotherms, and the size of the error bar in temperature is half of the diherence 
in temperature of these two isotherms. Finally, the density associated with the CMT is the 
critical micelle concentration at the CMT. A conservative estimate of the CMC at the CMT 
was obtained from the highest temperature simulated isotherm which contained micelles. 
Using this isotherm, the CMC at the CMT is within the density range between the CMC 
and high density boundary of the micellar fluid just below the CMT. 

Structural transitions at low temperature were determined using parallel tempering sim¬ 
ulations. Conventional WL-TMMC simulations in the grand canonical ensemble at hxed 
T roughly identihed the temperature region where elongated clusters formed, and where to 
perform the parallel tempering simulations. In the parallel tempering simulations, many 
isotherms from grand canonical WL-TMMC simulations were performed at closely spaced 
intervals in temperature to target the micelle-to-elongated cluster transition region, and the 
isotherms were allowed to exchange conhgurations between temperatures at constant num¬ 
ber of trimers. The conhguration exchanges in parallel tempering allowed for more frequent 
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sampling of the micelle to elongated cluster transition. This transition was identihed by 
both structural and thermodynamic dehnitions. In the structural dehnition, the transition 
occurs at the temperature at which there is equal probability of observing more than one 
micelle, and one elongated cluster. The maximum number of trimers, N^ax, in each isotherm 
was set to a value which would typically contain two micelles, when above the transition 
temperature. The thermodynamic dehnition is the temperature at which there is a peak in 
the constant volume heat capacity.— The constant volume heat capacity, Cy is computed 
as 


{Cv)nvt 


{U‘^)nvt - {U)nvt 
ksT^ 


( 6 ) 


where U is the potential energy and {...)nvt is a canonical ensemble average. The grand 
canonical ensemble average is then obtained from Equation HI In these parallel tem¬ 
pering simulations, twelve isotherms were simulated every = 0.25, in the range 

G [4.5,7.25]. In order to study the density dependence of this transition temper¬ 
ature, two sets of parallel tempering simulations were performed with different volumes, 
I//ct 3 = 729 and 5832. 


IV. RESULTS AND DISCUSSION 

We studied the phase behavior of a trimer huid as a function of bond length L and 
bond angle 0. The interesting Ending is that dramatic changes in phase behavior can be 
caused by simple changes in the geometry of the trimer. In particular, the phase behavior 
changes dramatically from macroscopic huid phase separation without self-assembly at low 
L to self-assembly without huid phase separation as L increases up to L = a. In the 
special in-between case of L = OAa, both huid phase separation and self-assembly occurred 
simultaneously, where the latter resulted in the formation of elongated clusters. 

A variety of self-assembled structures were observed for Lja = 0.4,0.75,1, as shown in 
Figure [2l In particular, two predominant types of self-assembled structures formed in the 
density range of interest in this study. The hrst type of self-assembled structures can be 
described as micelle-like spherical clusters. These micelles were of variable size, depending 
on both the state conditions and the shape parameters of the trimer model, shown in Figures 
[Si,[2)d, and The second type of self-assembled structure can be described as elongated 
clusters, shown in Figures [2li, |2^, and[2]F. One important feature of elongated clusters is that 
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FIG. 3. Fluid phase coexistence for (squares) L = 0, (circles) L = 0.25a, (triangles) L = O.Aa 
with 0 = 7r/3. The critical points, shown by the symbols at the maximum temperature, were 
obtained from fit to the law of rectilinear diameters. The labeled black circle corresponds with 
the structure shown in Figure [2]F. The error bars, smaller than symbols, were obtained from three 
independent simulations. 


they may form at low density. Also note that these elongated clusters may vary in shape, 
depending on both L and 0. A third type of self-assembled structure, packed cylinders, 
shown in Figures [2^ and [2)i, is only observed at high density. Simulations in this high 
density regime are beyond the scope of this study due to sampling difficulties. 

Fluid phase separation was observed for Lja = 0,0.25,0.4 and 0 = vr/S, as shown in 
Figure [3l In Figure [3], the density, p, is normalized by the excluded volume. Vex (see Table 
[I]) , in order for the phase coexistence curves for different values of L to be in a similar density 
range. For the largest bond length exhibiting fluid phase separation, L = 0.4a, both the 
low- and high-density coexisting liquids are inhomogeneous due to the presence of elongated, 
self-assembled structures (see Figures [2fe andEf). Incidentally, we observe an approximately 
linear dependence of the critical temperature on the bond length in the range investigated, 
as shown in Figure 01 The line in Figured] is a linear £t, fcBT/**(L)/e = aLja -\- b, where 
a = —1.40(6) and b = 0.86(2). This linear trend is reminiscent of a linear trend in the 
critical temperature with respect to relative attraction between two beads of a dimer model 
reported previously.— 

In contrast to the macroscopic phase separation observed for bond lengths L < 0.4a, 


12 






Trimer Self-Assembly and Fluid Phase Behavior 



FIG. 4. Critical temperature as a function of bond length with 0 = vr/S, shown by symbols, with 
linear fit. The error bars, smaller than symbols, were obtained from three independent simulations. 


trimer fluids with L > 0.75a self-assembled into micelles and did not exhibit macroscopic 
phase separation. In Figure [5l we show phase diagrams using the approach described in 
Section IIIIBI and Appendix for the following four trimer model parameter pairs denoted 
as (L, 0): (cr, vr/S), (cr, 7r/2), (cr, vr/d), (0.75 ct, vr/S). For all cases, the critical micelle concen¬ 
tration increased with temperature. In addition, the critical micelle concentration increased 
with bond angle, 0, for L = a at hxed T. At a select temperature, error bars for the 
CMC and the high density boundaries of the micellar fluid were computed as the standard 
deviation from three independent simulations with volumes, V/a^ = 512,729,857.375. A 
secondary purpose for using different values of V in the three independent simulations was to 
verify that the results were not system-size dependent, and indeed they were not dependent 
on system size. The critical micelle temperature increased with decreasing bond angle, 0, 
for L = a. For the spherical to elongated cluster transition temperature, a set of parallel 
tempering simulations yielded a small density range for the transition temperature. Recall 
that this represents a low-temperature boundary for the spherical micellar fluid. In order 
to investigate the density dependence of this spherical to elongated cluster transition over a 
greater density range, two sets of parallel tempering simulations were performed at volumes 
Vja'^ = 729 or 5832. The spherical to elongated cluster transition appeared to be relatively 
insensitive to density within the error bars of the simulations. 

In Figure [6|, the average cluster size as a function of temperature and density is shown 
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FIG. 5. Phase diagrams for (top left) L = a,Q = tt/S (top right) L = a,Q = 7r/2 (bottom left) 
L = fj, 0 = 7r/4 (bottom right) L = 0.75(7,0 = 7r/3. The labeled black circles correspond with 
the structures shown in Figure [2j The critical micelle concentration is shown with the (green x) 
structural definition and (black -|-) thermodynamic definition. The critical micelle temperature is 
shown by the red triangle. The micelle to elongated cluster transition is shown with the (open blue 
triangle) structural definition and (open blue square) thermodynamic definition. The high density 
boundary of the micellar fluid is shown by the solid red squares. Lines are guides to the eye. 

for the following four trimer model parameter pairs denoted as (L, 0): (a, tt/S), (a, 7r/2), 
((7, 7r/4), (0.75(7, 7r/3). Cluster sizes increase with increasing density and decrease with in¬ 
creasing temperature. In addition, cluster sizes decrease as the bond angle, 0, is increased 
from 7r/4 to 7r/2. Only temperatures above the spherical to elongated cluster transition are 
shown in Figure [6l This T range was chosen because, when elongated clusters form, the 
majority of trimers in the system are part of a single cluster, and the cluster size is trivially 
equivalent to the number of trimers in the system. 

The fluid phase behavior of the family of trimer models may be understood in terms of 
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FIG. 6. Average number of trimers in a cluster for (top left) L = a,Q = vr/S (top right) L = 
cr, 0 = Tr/2 (bottom left) L = a, 0 = 7r/4 (bottom right) L = 0.75cr, 0 = vr/S. Contour spacing is 
one trimer, and was obtained by 3-dimensional fit to a number of isotherms. The thick, solid black 
lines are the boundaries shown in Figure [5j The contours were truncated at 40 trimers. 

the relative size of the repulsive and attractive regions.— As the repulsive region and 
anisotropy increases, self-assembly is more favored than macroscopic phase separation. One 
way to quantify the relative sizes of the repulsive and attractive regions is to assume that the 
attractive region stays hxed while changes in the repulsive region are due to changes in the 
net excluded volume of the beads. It thus follows that the smallest bond length corresponds 
to the smallest repulsive region, and the longest bond length corresponds to the largest 
repulsive region. The models in Figure [T] and Table [I] are listed in order of decreasing ex¬ 
cluded volume. As the excluded volume increases with L, the trimer shifts from fluid phase 
separation at low L to self-assembly at high L, with the special case of L/cr = 0.4 possessing 
both fluid phase separation and self-assembly. The relative change in the size of attractive 
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to repulsive regions may also be observed via the Boyle temperature, /cs^Boyie/e- The Boyle 
temperature is the temperature at which the second virial coefficient is zero (see Section 
HID, where the attractive and repulsive contributions cancel to yield the compressibility of 
an ideal trimer fluid. Finally, the relative size of the attractive region with respect to the 
repulsive region also explains the change in the critical temperature, ksTc/e, where ks is 
the Boltzmann constant. Given that the trivial, isotropic case of L = 0 possesses a critical 
temperature, increasing the repulsive region reduces the critical temperature. Increasing 
the anisotropy of the model yielded self-assembled structures, but also reduced critical tem¬ 
peratures. Only in special cases may one observe both self-assembled structures and phase 
coexistence (e.g. Figures |2] and [3] and Refs. I88l-l36h. 


The results in this study are consistent with the experimental and computational results 
of Wolters et. ai, who studied a trimer fluid similar to the ones studied in this work.— The 
results in this study may also be used to guide future experiments on tuning the shape of the 
trimers to control the formation of self-assembled structures that have not been observed 
in experiments. While the attractive interaction in Wolters et. al. is shorter-ranged than 
in this work, the fluid phase behavior of their model can be cast within the context of 
this work by matching the second virial coefficient, which assumes Noro-Frenkel extended 
corresponding-states applies.—*^ Although we have not computed the phase diagram for 
the fluid studied by Wolters et. ai, we can anticipate its phase behavior by comparing the 
second virial coefficient, B 2 , with experiment, and also the Boyle temperature with Table 
m This comparison is made based upon a trimer model, referred to as the MM-LJ model 
in this work, using the same shape parameters of Wolters et. al, L = 0.57a, 0 = 91 
degrees and smaller repulsive bead sizes, ar = 0.85cr, but with the potential described in 
Section HI] of this work.— Wolters et. al. reported a second virial coefficient of B 2 /a^ ~ 
which corresponds to a deplentant concentration of cpd ~ 0.2. This value of B 2 /a^ in turn 
corresponds to a reduced temperature of ksT/e = 0.355 ±0.005 for the MM-LJ model. The 
results in this study are consistent with the possibility that the MM-LJ model forms only 
elongated clusters, as found in the work of Wolters et. ai, because the Boyle temperature 
for the MM-LJ model is fc^TBoyie/e = 0.735 ±0.005. This value of the Boyle temperature lies 
between models in Table |T] that both formed elongated clusters, but the models transitioned 
from forming spherical clusters at L = 0.75a to not forming spherical clusters at L = 0.4a. 
One possible conjecture is that tuning the trimer shape in experiments for increased repulsion 
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with respect to attraction (e.g. increased size of repulsive beads, ar, and L) may lead to the 
formation of spherical clusters. 


V. CONCLUSION 


The phase diagrams of trimer particles with one central attractive bead and two repulsive 
beads were computationally mapped out as a function of the trimer shape. It has recently 
been shown that it is possible to synthesize similar colloidal trimer particles,—*^ and this 
computational study may guide future experimental studies of different trimer geometries. 
The trimer particles self-assembled into spherical clusters, elongated clusters and packed 
cylinders. The shape of the trimers, and the state conditions, played a role in determining 
the type of self-assembled structure that is formed. In addition, some trimer geometries 
led to macroscopic fluid phase separation. The transition from microscopic self-assembly 
to macroscopic fluid phase separation may be understood in terms of the relative size of 
repulsion and attraction in the particle. In special cases, both self-assembled structures and 
macroscopic phase separation occurred simultaneously. 

While the effect of the shape of the trimer on the phase behavior is the emphasis of 
this study, future investigations may utilize interaction potentials that model a particular 
system more closely (e.g. shorter-ranged interactions for patchy colloids). Note that the 
continuous potential in this work was chosen to make the model amenable to molecular 
dynamics simulations, which will be the subject of future publications to study the kinetics 
of assembly. 


The obvious case of L = 0.5a was omitted from this study for the following reasons. 
Although extensive simulations were conducted for L = 0.5a, the state conditions where 
the fluid potentially exhibited phase separation and/or self-assembly involved low tempera¬ 
tures. Our existing set of Monte Carlo moves were not sufficient to sample these conditions 
adequately. Proper sanipling under these conditions may require more sophisticated cluster 
trial moves (e.g. Refs. I39I- I41I1 . Additional study of the cases near Lja = 0.4 and 0.5 may 
be the subject of future publications. 
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Appendix A: Examples of Determining Self-Assembly Structural Transitions 

In this appendix, examples of the structural transitions from the spherical micellar fluid 
are provided for select state points for L = cr, 0 = vr/S. As described in Section UlI B1 these 
transitions include the critical micelle concentration (CMC), the critical micelle temperature 
(CMT), the spherical micelle to elongated cluster transition, and the high density boundary 
of the micellar fluid. The transition between a free trimer fluid and micellar fluid is dehned 
by the critical micelle concentration (CMC). The CMC was obtained both structurally and 
thermodynamically, as described in Section UlIBl Figures [7] and E] illustrate these different 
approaches. In addition, there is a critical micelle temperature (CMT), above which a 
trimer fluid exists without micelles. This CMT is not a true critical point, but it is a 
useful construct that suffers from some arbitrariness. As demonstrated in Figure [9l micelles 
formed at fc^T/e = 0.275 due to the presence of a system-size dependent density of a second 
peak in the macrostate distribution.— But they did not form at the higher temperature of 
ksT/e = 0.3. Therefore, the CMT, Tcm, is between these two temperatures, and reported as 
ksTcm/e = 0 .2875±0.0125 in FigureO There is also a high density boundary for the micellar 
fluid, where micelles deform to improve packing, and eventually form different structures 
(e.g. cylinders shown in Figures |2^ and [2)i). This high density boundary was dehned 
approximately as the density at which the concentration of free trimers and premicellar 
aggregates is no longer constant. Finally, at lower temperatures, there is a transition between 
spherical micelles and elongated clusters. The structural dehnition of this transition was 
the temperature at which the system had an equal probability to form two micelles, or 
one elongated cluster. This is shown in Figure [10] as two spherical clusters of 20 trimers 
in size combined at lower temperatures. The thermodynamic dehnition of the micelle to 
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FIG. 7. Number density of free trimers and premicellar aggregates, Pfree with (black solid line) 
ksT/e = 0.25 and (red dashed line) ksT/e = 0.3. The bine shaded region shows where pfree is 
within 75 % of its maximnm value. L = a, & = vr/S, and V = 729it^. 



FIG. 8. Pressure as a function of trimer density for ksT/e = 0.25 with the ideal pressnre shown 
by the black dashed line and the fit to the second linear regime shown by the red dotted line. 
L = a, Q = 7r/3, and V = 4096(T^. 


elongated cluster transition was the temperature at which there was a peak in the constant 
volume heat capacity, shown in Figure [TT] For all cases, the thermodynamic and structural 
dehnition agrees within the error bars. The determination of the value of the spherical to 
elongated cluster transition temperature, and the errorbars, from two bracketing isotherms, 
are analogous to the precedure for the CMT described above. 
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FIG. 9. The probability to observe a number of trimers, 11 for L = a, Q = tt/S with (black 
x) ksT/e = 0.275, V = 729a^ (red circle) ksT/e = 0.275, V = 512a^ (black -|-) ksT/e = 0.3, 
V = 729(7^ and (red square) ksT/e = 0.3, V = 512cj^. When ksT/e = 0.275, n/ksT = 3ln{aA) — 
5. Otherwise, /i/Zc^T = 3ln{aA) — 4.5. A is the thermal de Broglie wavelength. Probability 
distributions are shifted by a constant for clarity. 



FIG. 10. Probability distribution of number of trimers in a cluster in parallel tempering simula¬ 
tions for L = fj, 0 = tt/3, V = 729(T^, and ksT/t = 0.154,0.16,0.167. 


Appendix B: Excluded Volume 

Excluded volume, Vex/<y^, was computed assuming hard sphere diameters of size a. In 
practice, the excluded volume was computed numerically by overlaying the trimer with a 
cubic grid of Np = 10® points and a side length, equal to a plus the maximum intra- 
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keT/e 


FIG. 11. Grand canonical ensemble averaged constant volume heat capacity, Cy, in parallel 
tempering simulations for V = 729cj^, L = a, Q = 7r/3, and N=35,36,37,38,39. 


particle distance from a site to the center-of-mass. By counting the number of grid points, 
Uo which overlap with the trimer, Vex = . By computing the excluded volume of a 

unit sphere and comparing to 47r/3, the numerical error is expected to be on the order of 
IQ-h 

Appendix C: Second Virial Coefficient and the Boyle Temperatnre 

The second virial coefficient, i? 2 (T), was calculated by Monte Carlo integration. 

B^(T) = -0.5 j^drf(r) = f(n) (Cl) 

/(r) = - 1 (C2) 

where Uj is the relative position and orientation of a second trimer with respect to the 
first trimer, and i = 1, randomly chosen positions and orientations of a second 

trimer with respect to the first. In practice, the cubic volume, V was chosen such that 
is greater than the twice the potential cut-off plus four times the maximum intra-particle 
distance from a site to the center-of-mass. Convergence was reached when |i? 2 |/c’'biock < 10“^ 
or cTbiock < 10“^, where (Tbiock is the standard deviation obtained from block averages of size 
Ntriai = 10®. The Boyle temperature, Teoyie, was found by starting at ksT/e = 0.15 and 
incrementing T by 0.01. The reported fcnTBoyie/^ in Table |T] was the average of the two 
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temperature increments nearest B 2 = 0, and the reported error was ±0.005 to span the 
entire temperature increment. 
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